Dissipative dynamics of atomic Hubbard models coupled to a phonon bath: 
Dark state cooling of atoms within a Bloch band of an optical lattice 
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We analyse a laser assisted sympathetic cooling scheme for atoms within the lowest Bloch band of 
an optical lattice. This scheme borrows ideas from sub-recoil laser cooling, implementing them in a 
new context in which the atoms in the lattice are coupled to a BEC reservoir. In this scheme, excita- 
tion of atoms between Bloch bands replaces the internal structure of atoms in normal laser cooling, 
and spontaneous emission of photons is replaced by creation of excitations in the BEC reservoir. 
We analyse the cooling process for many bosons and fermions, and obtain possible temperatures 
corresponding to a small fraction of the Bloch band width within our model. This system can be 
seen as a novel realisation of a many-body open quantum system. 

PACS numbers: 03.75.Lm, 32.80.Pj, 42.50.-p 



I. INTRODUCTION 



New frontiers in atomic physics have often been en- 
abled by the development of new cooling techniques. Ex- 
amples are provided by laser cooling and evaporative 
cooling [H, H H Ej|, which underly the exciting exper- 
imental advances to realize Bose Einstein condensates 
(BEC) and quantum degenerate Fermi gases of atoms 
and molecules @, i, 0, S H El EH El El El El ■ Re- 
cently, quantum degenerate gases of bosons and fermions 
have been loaded into optical lattices in one, two and 
three-dimensional configurations [H, EH El El HI HH 
I22I I23I HI [2f| . This makes it possible to realize atomic 
Hubbard dynamics with controllable parameters, open- 
ing the door to the study (and simulation) of strongly 
correlated systems with cold atoms. Control via external 
fields allows the engineering of atomic lattice Hamilto- 
nians for boson, fermion and spin models [H [2?], HI], 
which for a long time have been the focus of research in 
theoretical condensed matter physics. However, one of 
the most challenging obstacles in the realization of some 
of the most interesting condensed matter systems in cur- 
rent experiments is the need for lower temperatures of 
the atoms within a Bloch band of the optical lattice [53| . 

In this paper we will analyze a configuration where 
atoms moving in the lowest Bloch band of an optical lat- 
tice are cooled via laser assisted sympathetic cooling with 
a heatbath represented by a BEC of atoms. A unique 
feature of the present scheme is that the achievable tem- 
peratures of the atoms within a single Bloch band are 
significantly lower than those of the cooling reservoir, 
reaching a temperature of a small fraction of the Bloch 
band width. From a physics point of view, a guiding idea 
of the present work is formal analogies with laser cool- 
ing P, 0, H 0| ■ In laser cooling the motional degrees of 
an atom are coupled via laser excitation of the electrons 
to the effective zero temperature (photon) reservoir, i.e., 
the vacuum modes of the optical light field. The sponta- 



neous emission of photons carries away the entropy allow- 
ing a purification (cooling) of the motional state of the 
atom. In the present work we will follow a similar sce- 
nario by coupling atoms moving in the lowest Bloch band 
via laser assisted processes to an excited band, which can 
decay back to the lowest band by spontaneous emission 
of phonon (or particle-like) excitations in the BEC. We 
note that this cooling process is atom number conserving, 
in contrast to filtering and evaporative cooling tcchnquics 
[HIM HI- The present paper focuses on a specific proto- 
col for cooling within a Bloch band which follows analo- 
gies with (subrecoil) Raman laser cooling, as developed 
in seminal work by Kasevich and Chu [31| ■ 

We see an important aspect of the present work in 
pointing out the formal analogies between open quantum 
systems familiar from quantum optics with light fields 
and the present setup of atoms in optical lattices coupled 
to a phonon bath. These analogies not only stimulate the 
transfer of well established ideas of, e.g., laser cooling, to 
a new context, but also provides a new realization of an 
open quantum optical system with dissipative dynamics 
of cold atoms in optical lattices. The present paper ex- 
tends and gives details of the analysis of work published 
in Ref. (32| , in which setup we assume a homogeneous op- 
tical lattice, without an additional trapping potential for 
lattice atoms. As an additional remark, we also present 
an extension of these ideas to a form of spatial sideband 
cooling in a harmonic trapping potential. 
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FIG. 1: Setup for the Raman cooling scheme in the optical 
lattice, a) A Raman laser setup couples atoms in the lowest 
Bloch band to the first excited band, b) The lattice is im- 
mersed into a BEC of species 6, which serves as a cold reservoir 
for the lattice atoms, leading to decay of lattice atoms in an 
excited motional state back to the lowest band via collisional 
interaction with the reservoir and the creation of Bogoliubov 
excitations in the BEC (as sketched in the far right of the 
figure). 



as a reservoir for the lattice atoms, and should be seen 
as the counterpart of the vacuum modes of the radia- 
tion field in laser cooling, carrying away motional energy 
from the system. In our scheme, cooling is achieved in 
two steps: (i) We design a sequence of excitation pulses 
which efficiently excites atoms with high quasi-momenta 
in the lattice, but do not couple to atoms with q = 0. 
(ii) In a second step, the coupling to the BEC reservoir 
leads to decay of excited lattice atoms back to the ground 
state via emission of a Bogoliubov excitation into the 
BEC, thereby randomizing the quasi-momcntum. Re- 
peating these two processes leads to an accumulation of 
the atoms in a narrow region around q = 0, i.e., to cool- 
ing, in analogy to the Kasevich-Chu scheme of Raman 
laser cooling for free atoms [HI, [34|. We will show below, 
that our method works away from the limit of unit filling 
of the lattice, and is capable of cooling single atoms and 
many non- interacting bosons and fermions. The method 
can be utilized in a scheme to create strongly interact- 
ing phases of atoms in the optical lattice, by first cooling 
non-interacting atoms (thereby exploiting the tunability 
of interactions, e.g., via Feshbach resonances) and subse- 
quently ramping up the interaction. 

The paper is organized as follows. In the following sec- 
tion[TT]we will present a short overview of laser cooling as 
is relevant for understanding the physical analogies with 
our cooling method outlined in later sections. Sec. IIIII 
presents our model. We will then illustrate the details 
of the cooling protocol for the case of a single particle in 
Sec. IIVI In Sec. [V] we will discuss how the system can 
be adapted to achieve cooling for many non-interacting 
bosons and fermions and we investigate the affects of in- 
teraction on the cooling scheme for many bosons. In 
Sec. I VII we will show how the interaction strength be- 
tween already cooled lattice atoms can be ramped up to 
achieve cold strongly correlated gases. In Sec. lVIII we dis- 
cuss how similar ideas could be used to compact Fermions 
in a harmonic trap, and we then conclude and summarize 
the ideas and main results in Sec. IVIII1 



II. LASER COOLING 

In this section wc will outline the main ingredients 
needed for the description of laser cooling and give a 
short overview of the development of the different laser 
cooling schemes. The structure presented here is directly 
related to the cooling scheme we present in Sec. lIIIi as are 
the concepts of subrecoil laser cooling, which we discuss 
in more detail at the end of this section. 

As described in the introduction, in a laser cooling 
scheme, an atom interacts with a laser and is coupled 
to the electromagnetic radiation field, i.e., one considers 
a system of the form 



H = H a + Hb + Hj 



hit 



H 



LS- 



(1) 



Here, the atomic Hamiltonian H a , includes the atom's 
motion, possibly an external trapping potential and the 
internal structure, usually modelled as a two or a three 
level system. The atom interacts with the electromag- 
netic radiation field, described by Hi,, and this interaction 
is described by H lrA . In addition, a laser setup couples 
the atom's internal states, which is described by i?LS- 

The two steps of laser cooling, i.e., the upconversion 
of the energy due to the laser and the removal of system 
energy due to spontaneous emission, can be described in 
terms of a master equation in the Born-Markov approx- 
imation for the reduced system density operator p (see. 
e.g., |H), 

r 

p = i[H a + H LS ,p\- -(a + a p + pd+o ) 
r +i 

+ T J duN(u)e ik,iu &- pa+c- ik ' iu , (2) 

where we have specialized to atomic motion in one di- 
mension and to a two level system for sake of simplicity. 
The atom's position operator is denoted x, a + (<x~) is 
the raising (lowering) operator corresponding to the elec- 
tronic transition, T denotes the linewidth of the excited 
state, and ki the wavenumber of the laser light. Note 
that photons can be spontaneously emitted in all three 
dimensions, and the angular dependence of the sponta- 
neous emission is accounted for by the normalized dipole 
distribution N(u) in this one dimensional model. 

These basic ingredients have been utilized in the de- 
velopment of various different laser cooling schemes over 
the last years. This began with Doppler cooling 
where two counter-propagating laser beams incident on 
an atom, so that the radiation pressure of the two beams 
compensates for an atom at rest, but is Doppler shifted 
towards resonance and thus enhanced for the counter- 
propagating beam in the case of a moving atom. In 
Doppler cooling the final temperature is limited by the 
linewidth of the excited state. Lower temperatures were 
obtained in various schemes, including polarisation gradi- 
ent cooling and Sisyphus cooling [361 ] , where the limiting 
temperature is given by the recoil energy an atom re- 
ceives during the emission of a single photon. Cooling 
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of atoms even below this single photon recoil limit was 
proposed and observed in the form of Raman laser cool- 
ing and velocity selective coherent population trapping 



A. Subrecoil Laser Cooling 

The basic idea in subrecoil laser cooling is to make 
the photon absorption rate velocity dependent and, in 
particular, vanishing for a dark state. The atoms will 
consequently accumulate in this state during the cooling 
process. In the following wc will describe the examples of 
Raman cooling and velocity selective coherent population 
trapping (VSCPT). 



Cooling is achieved in a two step process. In the first 
step a set of Raman laser pulses is designed that effi- 
ciently transfer atoms with high momenta in the ground 
state <7i to (72, but do not couple atoms with zero mo- 
mentum. In the second step, which is applied after each 
single excitation pulse, atoms are optically pumped back 
to the state g\ via the excited state e (see Fig. [5]) . The 
spontaneous emission process randomizes the momentum 
of the atom, which leads to a finite probability of falling 
into a region near the dark state, which in this case is the 
state I g\ , q = 0) . Repeating the two steps leads to accu- 
mulation of atoms in state | g\) with momentum near 
zero, i.e., to cooling of the atom. 

2. Velocity Selective Coherent Population Trapping 
(VSCPT) 
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FIG. 2: Schematic picture of one excitation and decay step in 
subrecoil Raman laser cooling, a) Atoms from a region with 
high momentum \q\ > are transferred from one ground state 
gi to a second ground state <?2- b) The atoms in state gi are 
optically pumped back to the initial state g\ via the electron- 
ically excited state e, thereby randomizing the momentum 
due to the recoil kick. Cooling is achieved by designing a 
sequence of excitation pulses to efficiently transfer all atoms, 
except for those in the dark state with q = 0, to the second 
internal state, each pulse is followed by the optical pumping 
process, and repeating the sequence leads to accumulation of 
atoms in the dark state, i.e., to cooling. 



In Raman cooling [3l| one considers a A-system as 
shown in Fig. [5J consisting of two (hypcrfme) ground 
states g\ and g^ and one electronically excited state e. 
The energy of the moving atom is given by the sum 
of its internal energy u>i, where i = gi,g2,e, and its 
kinetic energy q 2 /2m. The state \gi,q) is coupled to 
I g-2, q+ (k\ + k-2)) via a Raman process Z/ls, which is 
far detuned from the excited state. Here, q denotes the 
momentum of the atom, k\ and ki are the wave numbers 
of the two Raman laser beams and we use units where 
h = 1 throughout the paper. 




FIG. 3: a) Schematic picture of the setup in VSCPT. Two 
circularly polarized laser beams incident on the atom, couple 
the two ground states \g~i) and \gi) via the excited state 
I eo). A dark state (see text) appears due to quantum inter- 
ference, where atoms accumulate during the cooling process 
as they decay there via spontaneous emissions from the ex- 
cited state, b) Theoretical model for the excitation rate in 
the Levy statistics analysis. For momenta \q\ < qo the excita- 
tion rate follows a power law, outside this region it is assumed 
constant. 



In VSCPT (see Fig.[3Ji)) [3J|, two counter-propagating 
circularly polarized laser beams couple the two ground 
states and | gi) via the excited state |eo). In the 

setup states with angular momentum J = 1 are used 
in the ground and excited state manifold, so that only 
the three states | g-i), \ gi) and | eo) are coupled via the 
lasers, and spontaneous emission only leads to decay back 
to the states g-\ and g\ since the transition eo — > go is 
forbidden. A dark state forms due to quantum interfer- 
ence and is given by | D) = {\g^,ki) + | gi, -k{))/y/2, 
where ki is the wave number of the two lasers. This dark 
state does not couple to the Raman process and will be 
increasingly populated during the cooling process due to 
the decay of atoms in the excited state via spontaneous 
emission. This leads to cooling of the atom, characterized 
by two narrow peaks at k\ and —ki in the final momen- 
tum distribution of the atom, as the dark state | D) is a 
superposition of these two momentum eigenstatcs. 
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B. 



Levy statistics 



tribution and find 



Analytical calculations based on Levy statistics [37j | 
have been shown to be a very powerful and accurate tool 
in the context of subrecoil laser cooling. They will be ap- 
plied below to describe our optical lattice Raman cooling 
process. In this section we briefly review the underlying 
model and some results in the context of Raman laser 
cooling for the most important physical parameters, such 
as the temperature and the fraction of the momentum 
distribution in the dark state. 

In Levy statistics [37| a trapping region near q = 0, 
with | q | < gtrap and a recycling region with \q\ > <7trap is 
defined. The parameter gtrap is an auxiliary variable and 
the results for real physical quantities do not depend on 
it. During the cooling process, the atom will undergo a 
random walk in momentum space, with A/" trapping peri- 
ods of duration r, , alternating with J\f recycling periods 
of duration ft , during which the atom is in side and out- 
side the trapping region respectively. The total time 
of the cooling process can be written as the sum of the 
total trapping time TV and the total recycling time TV 

© = T^ + fv, TV = ^2 n, TV = ^2n, (3) 

i=l i=l 

where the r, and the Ti are independent random variables. 

The efficiency of the cooling process can be quanti- 
fied by the statistics of the total trapping and recycling 
times, which are determined by the excitation rate R(q) 
of atoms with momentum q from the state g\ to the state 
g-2- The rate R(q) is modelled as 



(6) 



R(q) 



1 

TO 



M < Qo, 



R(q) = -, \q\>qo, 

TO 



(4) 
(5) 



as schematically plotted in Fig. [3)d). The values of ro, 
<?o and A depend on the details of the excitation pulses 
and are mainly determined by the duration, momentum 
transfer and shape of the last two Raman pulses. In 
subrecoil cooling, typically A > I (especially A = 2 for 
the case of time square excitation pulses in Raman laser 
cooling) and thus the probability distribution for the to- 
tal trapping times TV is a broad distribution, i.e., the 
expectation values (TV) and (Tj^-) diverge. In this case, 
a generalized central limit theorem (for details see 
predicts a Levy distribution for the probability distribu- 
tion for the total trapping time. 

This distribution is the starting point for a calculation 
of the relevant physical parameters, like the width of the 
momentum distribution or the height of the central peak 
at q = (for details again see 37] ), with the following 
results. We define the temperature in terms of the half 
width Aq at e -1 / 2 of the maximum of the velocity dis- 



especially we find T cx O -1 for the case of time square 
excitation pulses. Similarly we derive an expression for 
the population density at q = and time 0, no(O), which 
is given by 

, . A 2 sin7i7A 1 /6\ 1/A m 



27T7(1/A) qo \to 



where in this equation, j(x) denotes the Euler gamma 
function. 



III. RAMAN COOLING IN AN OPTICAL 
LATTICE 

In this section wc will describe our Raman cooling 
scheme for atoms in an optical lattice, therefore we will 
again shortly describe the idea of the cooling process and 
then analyze the different parts of the setup in detail. 

In our scheme cooling of the atoms is achieved in a 
two steps (c.f. Figs. Hk) and HJs)) and can be seen in 
analogy to free space Raman cooling (Fig. [J), (i) Ra- 
man laser pulses are designed to excite atoms with large 
quasi- momenta \q\ > from the lowest band to the first 
excited band, whilst not exciting atoms in the dark state 
(with q — 0) (sec Fig. BJi)). (ii) The decay of excited 
lattice atoms goes along with the emission of a phonon 
(Bogoliubov excitation) into the BEC reservoir (see [4})). 
Assuming a BEC temperature /cbT, -C uj the Bogoliubov 
modes with an energy corresponding to the separation of 
the Bloch bands are essentially in the vacuum state and 
the BEC effectively acts as a T = reservoir. This is 
in analogy to the coupling of the atoms to the vacuum 
modes of the radiation field, giving rise to spontaneous 
emission into a T = reservoir in the case of laser cool- 
ing. Repeating these two processes leads to cooling of 
the atoms, as they accumulate in a region near the dark 
state with q = 0. 

Note, that in contrast to the free space version of Ra- 
man cooling, in our setup the upconversion of energy is 
already performed in the excitation step. This is only 
possible, because in our setup the spontaneous decay 
back to the lowest band is due to collisional interaction 
with the BEC reservoir which can be switched off dur- 
ing the excitation step, e.g., via Feshbach resonances. 
A direct analogue to the cooling protocol in free space 
Raman cooling could, however, also be achieved with a 
slightly different setup: For example, one could use a 
spin-dependent lattice, and perform the excitation step 
from the lowest Bloch band for lattice atoms a to the 
lowest Bloch band of a second species a', followed by an 
optical pumping step via an excited Bloch band, as in 
free space Raman cooling. 
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In our setup we assume a homogeneous optical lattice, 
i.e., no additional external (harmonic) trapping poten- 
tial for the lattice atoms. This assumption is an exper- 
imental requirement inherent in the realization of many 
strongly correlated phases, and a topic of significant cur- 
rent interest experimentally (see, e.g., (38^). Our analy- 
sis will be performed for the ideal case of a homogeneous 
system, and inhomogcneitics appearing in a real exper- 
iment will limit the achievable temperatures. In prac- 
tice, we expect that advances in homogeneous traps and 
cooling methods will occur somewhat iteratively, in that 
clean flat-bottomed traps allowing better cooling, which 
will provide the opportunity to cool atoms to even lower 
temperatures where the remaining imperfections become 
more noticeable. 




7r qd 



FIG. 4: Schematic picture of one excitation and decay step in 
Raman cooling in an optical lattice, a) Atoms are transferred 
from a region with high quasi-momentum \q\ > in the low- 
est Bloch band to the first excited band, b) The collisional 
interaction with the BEC atoms is switched on, the resulting 
decay of the excited lattice atoms leads to a randomization of 
the quasi-momentum. Sequences of pulses, each one followed 
by a decay time t c , efficiently excite all atoms outside a nar- 
row region around q = 0. Repeating the sequence leads to 
accumulation of atoms in the dark state region around q = 0, 
i.e., to cooling. 



A. Lattice atoms and laser setup 

The Hamiltonian H a for the lattice atoms can be writ- 
ten as H a = Hq + Hj , where 



H a = 



R 3 



d 3 x #( x ) + V a (x)) 4(x), 



(8) 



describes the kinetic energy of the atoms and the 
optical lattice potential V a (x) = V a>x sin 2 (7ra/d) + 
V at y sm (ny/d) + V a , z sm 2 (-Kz/d). Here, d = A/2 is the 
lattice spacing with A the wavelength of the lasers gener- 
ating the lattice potential, and V a j is the strength of the 
optical lattice potential in j — x, y, z direction, m a is the 
mass of the lattice atoms and V4(x) and V> a ( x ) are the 
field operators for the lattice atoms, which will satisfy 



(anti-)commutation relations in the case of (fermions) 
bosons. Onsite interactions between lattice atoms are 
represented by 



H, 



9a 



d 3 x ti(x)iPi(x)Mx)Mx), (9) 



2 Jn 3 

^Tra aa /m a with a aa the s-wave scattering 



where g aa 
length. 

The interaction of the two laser beams generating the 
Raman transition with the lattice atoms in the rotating 
wave approximation can be written as 



#LS = 



d 3 a 



(( 



-h.c. 



^i( x ) , gaW ist 
2 + 2 

A#( x )^ e ( x ) 



(10) 



where the Raman detuning S = uj\ — lo 2 is given by the 
frequency difference of the two lasers, A is the detuning 
from the excited state, and fii and £^2 denote the single 
photon Rabi frequencies of the two lasers. Both lasers are 
far detuned from resonance with the transition between 
ground and excited state, which can thus be adiabatically 
eliminated (see Appendix [A"|l . 

For simplicity, wc will assume in the following an 
anisotropic lattice potential, with V a<v = V Qj2 = V a ,i_ ^> 
V aiX so that the atoms effectively move in one dimension 
along the ^-direction, whereas the transverse hopping is 
suppressed due to the large trapping potential. Such an 
analysis is readily extended to higher dimensions. In the 
case where only the two lowest bands play a role in the 
dynamics of the system, we can express the field opera- 
tors in terms of the Bloch functions for the lowest two 
bands in the x-direction, <p^(x), where the band index 
a G {0, 1}, and localised wavefunctions Wy(y) and w®(z) 
representing the confinement in the transverse directions. 
This imposes requirements on the laser coupling strength 
and detuning, and also on the interaction strength be- 
tween atoms in the lattice (see below). We obtain 

to= E £<(*K(yV°(^K. (ii) 

a=0,l q 

where (^4^)^ creates a particle with quasi-momentum q 
in Bloch band a. The operators (A^y and (A^) will 
again satisfy (anti-)commutation relations for (fermions) 
bosons. Inserting this into Eqs. ©-([H]), we obtain 

Ho + ^ls = E e % ( k i) f K + - s ) E ( A l) f A l 



(/.a 

o 
2" 



Ef( i O ti °-* +h - c - 



(12) 



where w is the energy separation of the Bloch bands, and 
the kinetic energy of the lattice atoms 



-2J a cos(gd). 



(13) 
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In a setup where the wave number difference of the two 
running wave laser beams inducing the Raman process is 
parallel to the rr-direction and of magnitude 5q, we can 
define the effective Rabi frequency as 



n 



dcexp(— iSqx)w l (x)w° {x) . (14) 



Here, the Wannier functions w a (x) are defined as 



volume V by 



(19) 



where gbb = 47ra;,£,/m£,, with abb the scattering length for 
the interaction of reservoir atoms, and the operator 
creates a Bogoliubov excitation with momentum k and 
energy 



W a (x — Xj ) 



(15) 



where Xj is the position of the j-th lattice site with 
j G {1,..,M} and M is the number of lattice sites, 
the discrete quasi-momcntum q in the finite lattice qd € 
{— (M - l)/2, . . . (M - 1)/2}2tt/M. Note that we ne- 
glect contributions involving Wannier functions in differ- 
ent lattice sites. 

Similarly, we find that the terms describing onsite in- 
teractions between lattice atoms can be expressed as 



- E 

91 ,92, 93, a 



1+92-93 



M ^ 

91192,93 



u aa (A«y (Ai) 1 A« 3 Ai_ 
u 10 (Al)\Aiy A° 3 Al i+q2 _ qs , (16) 



with 



U aa =g ± l <k\w a {x)\ 2 \w a '{x)\ 2 . (17) 



Here, g± = 2uj±a aa and u>± = y/4:Va,±uiR is the oscilla- 
tion frequency of the transverse confinement, ojr is the 
recoil frequency. Note that this expression is valid for 
U aa <C u>. In order to obtain the expressions for eqs. (flU)) 
and (fT7)) we have inserted eq. (fTT|) into cq. ([9]) and per- 
formed the integrals over the transverse directions, where 
we have approximated the lattice potential by harmonic 
oscillators with frequency u>±. 

The interaction strengths (j 1 T[) can be explicitly calcu- 
lated if we also approximate the Wannier functions in 
^-direction with harmonic oscillator wave functions, and 
we find 



U 00 = g M / 1 ^(V a , x ui R y , 



(18) 



and U 10 = U 00 /2, U 11 = 3C/ 00 /4. 

In summary, the above model requires J a , U a ' a , il <C 



B. BEC reservoir and interaction with the lattice 
atoms 

The BEC-reservoir is described as a 3D homogeneous 
quantum gas consisting of Nb particles of mass nib in a 



Ik|4 



(20) 



in the reservoir. The sound velocity in the BEC is given 
by c = yj gbbPb/rrib, where pb = Nb/V is the mean con- 
densate density. 

The interaction of the lattice atoms with the supcrfluid 
reservoir is modelled by the density-density interaction 
Hamiltonian 

£int = 9ab f d 3 ^t(x)^(x)^(x)^ a (x), (21) 
JTR 3 

with the interaction strength g a b = 47ra f,/2m r , where 
a a b is the inter-species s-wave scattering length and m r = 
m a mb/ (jn a + TUb) the reduced mass. The field operators 
for the BEC can be expressed as 



(22) 



and 



Hi 



« = E( u ^ k c ikx + , k 6[e- ikx ) (23) 



in terms of creation and annihilation operators 6^ and 
6k for a Bogoliubov excitation with momentum k = 
(k,k y ,k z ). The coefficients Mk and Uk can be written 
as 



1 L k 

= -7===, W k = 



where 



Lk = (Ek - k 2 /2m b - m b c 2 )/m b c 2 



(24) 



(25) 



We neglect the terms proportional to Sipl{x.)Sijjb('x.) [33[ 
and using eq. (jllj) and leaving out the constant mean 
field shift we obtain 



H ^ = E E ( G ^k (A^i A a q '_ k + h. 



a, a' k,g 

Here, the coupling 



(26) 



^a,a' ~ 9al 



1/2 



dxe ikx w a (x)w a \x), 



(27) 
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where 5(k) is the static structure factor 



S(k) = (u u + v k ) 2 = 



2m h E v 



(28) 



In eq. (|2T|) we have again neglected the overlap of Wan- 
nier functions in different lattice sites and performed the 
integration over transverse lattice directions, where the 
wave functions are again approximated with ground state 
harmonic oscillator wave functions and result in factors 
of one, provided mbLv/2m a uJ± <C 1, which is usually ful- 
filled for our one dimensional lattice potential. The in- 
teraction Hamiltonian Eq. (f2"o]) describes scattering pro- 
cesses, where a Bogoliubov excitation with momentum 
k in the (three dimensional) BEC reservoir and a lattice 
atom with momentum q — k (with k the component of 
the momentum of the Bogoliubov excitation along the 
x-direction) in Bloch band a 1 are annihilated (created) 
and a lattice atom with momentum q in Bloch band a is 
created (annihilated). 

The structure factor S'(k) — > 1 for energies much larger 
than the chemical potential fi = m^c 2 = gbbPb, where 
the corresponding Bogoliubov excitations are particle- 
like, whereas for energies much less than \i the excitations 
are phonons and 5(k) oc |k| — > 0. In our setup we will 
choose 4J° </i<w (see Fig. [S]), which corresponds to 
typical experimental parameters. In this case, interband 
transitions will involve absorption and emission of Bogoli- 
ubov excitations in the particle branch (E^ ~ uj). These 
will have much larger coupling strengths (S'(k) ~ 1) 
than intraband transitions, which involve absorption and 
emission of Bogoliubov excitations in the phonon branch 
(E k ~ J Q , 5(k) oc |k| -> 0) (c.f. Fig. ED. 




FIG. 5: Overview over the energy scales in the system. Left 
part: Energy structure in the optical lattice in the quasi- 
momentum picture. Right part: The dispersion relation of 
the Bogoliubov excitations in the (three dimensional) BEC 
reservoir. Excitations corresponding to the band separation 
CO are particle-like, small excitations with energies on the or- 
der of a Bloch band width are phonon-like. The chemical 
potential is larger than the width 4J° of the lowest band, but 
less than the band separation to. 



In the derivation of the master equation (section fill C| 
we will restrict ourselves to interband processes, i.e. to 
the decay from the first excited Bloch band back to the 



lowest band. In the following we will comment on intra- 
band processes, which can lead to sympathetic heating 
or cooling of atoms within a Bloch band due to collisions 
with Bogoliubov excitations in the BEC reservoir. 

The heating/cooling process is described by the inter- 
action Hamiltonian ([2"o]) and the corresponding rate can 
be estimated, e.g., with Fermi's golden rule. The heat- 
ing/cooling processes involve the scattering of a lattice 
atom with a single Bogoliubov excitation, and Fermi's 
golden rule implies both energy and momentum conser- 
vation. We will illustrate this for the case of sympathetic 
heating within a Bloch band, as this represents a possi- 
ble imperfection for the cooling process (whereas sympa- 
thetic cooling would only speed up the cooling) and the 
same arguments are valid for intraband cooling. The typ- 
ical case of such a heating process in our scheme would 
be a scattering of an atom with momentum 5 « and 
energy £®~q to a momentum q' with energy s®, + c|k| 
(excitations within the Bloch band are sound waves), 
as typically most of the atoms are within a narrow re- 
gion around q = after a few iterations of the cooling 
process. This process is however only allowed if energy 
conservation c|k| = e° q — e q , and momentum conserva- 
tion along the lattice axis k = q — q' are fulfilled. As 

clkl = 



k^ + > c\k\ we find that this process 



is forbidden if the condition 



J°< 



y/fjiu R m a /(2m b ) 



(29) 



is fulfilled (see Fig. For the typical parameters in 
our setup (see Fig. 0, where (j,,u>r S> J , this condition 
will always be fulfilled, and thus sympathetic heating and 
cooling within a Bloch band are forbidden by energy and 
momentum conservation. Higher order processes arising 
from scattering with two or more BEC excitations will 
be small 1331. 



e(qd), E(k) 




tv qd,k 



FIG. 6: Schematic picture of the criterion for the occurrence 
of sympathetic heating/cooling within the lowest Bloch band. 
In analogy to the Landau criterion for superfluidity, the cor- 
responding scattering processes are only allowed if the Bloch 
band e(q) intersects the cone showing the Bogoliubov excita- 
tion energies corresponding to the relevant momenta k. 



8 



C. Master equation for decay 

The collisional interaction (f2"u| of atoms in the excited 
Bloch band of the optical lattice with the BEC reser- 
voir leads to a decay of excited lattice atoms back to the 
lowest band, in analogy with spontaneous emission (see 
Sec. [nj. F° r typical BEC temperatures fceTb <C uj, the 
Bogoliubov modes corresponding to the band separation 
u> will initially be in the vacuum state, and we can derive 
an effective T = master equation in the Born-Markov 
approximation for the reduced system density operator 
p, which describes atoms moving in the lattice, while 
the BEC is treated as a reservoir of Bogoliubov excita- 
tions. This is done in close analogy to [1^ and [33| in the 
context of lattice loading of Fermions from an external 
reservoir, and cooling of single atoms in a harmonic trap 
immersed in a BEC, respectively. We find (see Appendix 

u 



dr 2 



E 



pc[c k 



(30) 



where the one dimensional momentum k along the lattice 
axis is bounded by |fc| < /c max = y/2myjj due to energy 
conservation, and the jump operators are defined as 



Cfe 



5>S)t(aJ)e- 



ikx j 



q 



(31) 



with the position space operators af = 
(l/\/M) J2 q exp(igxi)A^. Note that an operator written 

Aq_ k should always be understood as A^,, where q' is 
a quasi-momentum in the first Brillouin zone, found by 
subtracting an integer multiple of the reciprocal lattice 
vector from q — k, i.e., q 1 = q — k — zG, z € %, G — 2ir/d. 
Similarly, the operator = Ck+ Z G- 

If a single atom is present in the excited Bloch band, 
the total decay rate via creation of excitations with all 
possible values of k is given by T = Tfc. We can then 
define the distribution of emitted excitations, dr/dfc, 
which for deep lattices (where we can approximate Wan- 
nier functions by harmonic oscillator ground states, and 
to ^> | J x |, J°) can be written explicitly as 



dT 

dF 



2tt 



gj b Pbm a a^k 2 a 2 fc 2 /2 



4tt 



(32) 



Here ao = ^/l/m a uj is the size of the ground state of the 
harmonic oscillator in the x-dircction and L = Md is the 
length of the lattice along the x-direction. 



The total decay rate can be explicitly calculated in the 
harmonic oscillator approximation as 



9abPb™-b 

27rao 



2^e- 




(33) 




^/2m, b u) 



3 d 2 d d d 2 d 3 3 



FIG. 7: The distribution of k values in decay events, dr/dfc, 
plotted in dimensionless units for the case where fc max ^> ir/d. 
In the plot we indicate integer multiples of the Brillouin zone 
width with dashed vertical lines. In order to compute the 
change in quasi-momentum of a lattice atom q' — q, the values 
for k must always be translated back into the first Brillouin 
zone. 



where erf(x) denotes the error function. The value of 
r can be tuned by changing the scattering length, the 
density of the BEC reservoir or the depth of the optical 
lattice. 

If we compute the distribution for the change in quasi- 
momentum, q' — q of the lattice atoms, we must translate 
values for k outside the first Brillouin zone, |fc| > n/d 
back into the first Brillouin zone, as q' = q — k — zG, z £ 
7L,G = 2iv/d. This is indicated by the regions between 
the dashed lines in Fig. [7] We can distinguish two inter- 
esting limits for our parameters based on the ratio of the 
upper bound on |fc|, k max = ^/2mjw and the extent of 
the first Brillouin zone, n/d = ^2m a ujR: 

(1) fcmax 3> n/d. Here, the distribution dT/dk ex- 
tends over k values much larger than the first Brillouin 
zone, as depicted in Fig. [7] When we compute the dis- 
tribution of changes in quasi-momentum for the lattice 
atoms, this will be approximately uniform over the first 
Brillouin zone. Note that this limit also corresponds to a 
situation in which the wavelength of emitted excitations, 
At = 27r/fc is typically much shorter than the separation 
between lattice sites d. Thus, collective effects of decay 
on different lattice sites (analogous to super-radiance and 
sub-radiance in atomic decay) are suppressed. 

(2) fcmax < it/d. Here, the distribution dr/dfc is cut off 
before it reaches the edge of the first Brillouin zone. As 
a result, the distribution of Aq is localised at low values, 
peaking at the cutoff value k = ±fc max . This can be used 
to target the decay to one area within the first Brillouin 
zone. For example, we can choose to excite atoms to 
a quasi-momentum value in the first band from which 
decay into the dark state will be strongly favoured. Note 
that this limit also corresponds to a situation in which 
the wavelength of emitted excitations, > d. Thus, 
collective effects involving decay on different lattice sites 
occur, and this decreased spatial resolution of the decay 
process corresponds to the increased resolution that we 
observe in momentum space. These effects are properly 
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accounted for in our calculations. 



IV. SINGLE PARTICLE COOLING 

In this section we will analyze the cooling process con- 
sisting of the two steps (i) the Raman laser excitation and 
(ii) the decay of excited lattice atoms. We will describe 
how efficient excitation laser pulses can be designed and 
present the results obtained from both numerical simu- 
lations and from analytical calculations based on Levy 
statistics. In the excitation step, which represents the 
first part of our cooling protocol, we will assume that the 
interaction with the BEC reservoir can be switched off 
(e.g., via a Feshbach resonance [iol. lill. [4^| ) . whereas the 
Raman coupling is switched off during the decay step. 



A. Designing the required laser pulses 





p(q) i 

0.5 


1 

0.5 


1 

0.5 



FIG. 8: A typical efficient pulse sequence consisting of three 
time square pulses, a) A short pulse excites atoms with large 
quasi-momentum at the edges of the Brillouin zone (dotted 
line), dashed line in b) and solid line in c), longer pulses excite 
atoms with quasi-momentum closer to zero, always keeping 
P(q = 0) = 0. Parameters used: Q = (27.4, 8.4, 8.4) J°, Sqd = 
(0.31, 2.12, -2.12), (5 - w) = -(28.2, 25.2, 25.3) J° 




We define the probability Pj{q) that the j-th pulse 
(with j = Q..N P — 1 and N p the number of pulses) excites 
an atom with initial quasi-momentum q from the lowest 
band to the excited band and require Pj (q) = for q f=a 
and Pj{q) — > 1 for states with higher quasi-momentum 
(c.f. Fig. 0Ji)). The probability Pj(q) can be obtained 
by solving the Hciscnbcrg equations of motion for the 
system operators, 



"l — — -n-,, 



(*). 



(34) 



where the effective detuning S q+ g q . = uo + £ \ + sq ~ £ q — 

These equations are valid in the subspace of a single 
lattice atom, where the interaction Hamiltonian Hi = 
0, and can be analytically solved for two simple cases. 
In the case of weak excitation, J Slj(t)dt < 1, we find 
probability 



in the case of a time square pulse, for the resonant tran- 
sition and adjust the parameters to always keep the first 
node of the sine- function at q = 0. Such a pulse sequence 
consisting of N p ~ 3 pulses, is shown in Fig. [8] We start 
with an intense laser pulse to excite atoms with momen- 
tum qd ~ 7T around the edges of the Brillouin zone, and 
then move the resonance closer to q = by adjusting the 
Raman detuning Sj and momentum kick Sqj and at the 
same time reducing the intensity of the laser beams. To 
be able to resolve the band structure with our Raman 
pulses we always have to fulfill SI <C 8IJ 1 ! and conse- 
quently t ^ 7r /8 1 1/ 1 1 . Note, that the relevant energy 
scale is the hopping | J x \ in the upper band, which is typ- 
ically an order of magnitude larger than the hopping J° 
in the lower band for lattice strengths V a , x ~ IQ^r- The 
parameters here have also been carefully chosen to avoid 
unwanted excitation to higher bands. 



B. Results 



PM = 



(35) 



in terms of the Fourier transform of the Raman Rabi 
frequency £lj(t). In the case of a time square pulse 
(£lj(t) = for < t < Tj and £lj(t) = otherwise), 
we have 



0$ 



Sill 



+ %/2. (36) 



The goal of the excitation step is to design effi- 
cient laser pulses, which excite atoms with large quasi- 
momentum |g| > but do not couple the atoms in the 
dark state q = 0. We would also like to do this on a fast 
timcscalc and therefore require a 7r-pulse, i.e., Tj = Tv/Qj 



In this section we will quantitatively analyze the cool- 
ing process, computing the final temperature and cooling 
timcscales. We make use of both numerical and analyti- 
cal methods, and compare the results we obtain in each 
case. 

In the numerical analysis we simulate the time evolu- 
tion of the system density operator using a monte carlo 
wave function method [43l . |44| . In the simulations we 
start with an initial mixed state according to a thermal 
distribution of atoms in the lowest Bloch band with a typ- 
ical temperature 4J° <C feT <C cj. To obtain good ap- 
proximations of the system density operator at all times 
we evolve the state according to the master equation (f3TJ)) 
and typically take the statistical average over ~ I0 5 tra- 
jectories in the simulations, for a one-dimensional optical 
lattice, with M = 101 lattice sites and two bands. 

The analytical calculations make use of Levy statistics, 
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similar to the corresponding calculations for freespace 
subrecoil laser cooling (see section Hi B|) . We define the 
temperature of the atom in terms of the width of the 
quasi-momentum distribution, and find 



k B T = 2J°{Aqd) 2 



(37) 



with Ag denoting the half width of the momentum dis- 
tribution at e -1 / 2 of the maximum. This is in close anal- 
ogy to free space Raman cooling (37| . with the mass in 
free space replaced with the effective mass in the optical 
lattice, and the quasi-momentum now playing the role 
of momentum in free space. We again find T oc _2//A 
(see section Hi B | and consequently time square excitation 
pulses again lead to efficient cooling, as A = 2 for time 
square pulses also in the presence of the optical lattice. 
Similarly we again find no(O) oc Q 1 ^. 

In Fig. [S] we compare the analytical and the numeri- 
cal results for the temperature and the fraction of atoms 
in the dark state. In the simulations, repeated applica- 
tion of the N p excitation and decay steps leads to the 
development of a sharp peak in the momentum distribu- 
tion n°(q) in the lowest Bloch band already after a few 
iterations, as can be seen in Fig. [9] a). 




5X10 1 tJ° 



FIG. 9: a) Development of the narrow momentum peak after 
0, 1, 3, 5 and 10 iterations of the pulse sequence shown in 
Fig. [5] for M = 101 lattice sites, b) Temperature in units of 
4J° against time in units of 1/J°. c) Fraction of the atoms in 
the dark state with q — against time in units of J°. In both 
plots, circles show numerical results obtained from quantum 
monte carlo wave function simulations, the solid line shows 
analytic results obtained from Levy statistics. Parameters as 
in Fig. [8] again for M = 101 lattice sites. 
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FIG. 10: a) Temperature and b) dark state population after 
a time t = 70 J for the parameters as in Fig. [5] as a function 
of the residual decay V obtained from numerical simulations 
as described in the text. 



temperatures of ksT/AJ ~2x 10 -3 and a typical frac- 
tion on the order of a few tens of percent in the central 
peak can be reached in tf J° ~ 50 for the pulse sequence 
shown in Fig. [8] In the simulations we have furthermore 
used r = Iujr, which can be obtained from eq. (|33p for 
Pb = 5 x 10 14 cm -3 and a s = 350a;, (e.g., via a Fesh- 
bach resonance) for 87 Rb in the lattice and 23 Na in the 
reservoir, a& is the Bohr radius. Note, that for the param- 
eters we use (see caption in Fig. [BJ), the cooling timescale 
is mainly determined by the duration of the last two ex- 
citation pulses. The cooling timescale is typically much 
faster than the excitation pulses, and thus smaller values 
of r docs not significantly slow down the cooling process, 
unless the cooling timescale becomes comparable to the 
duration of the excitation pulses. 

In our cooling scheme we assume that the interaction 
of lattice atoms with the reservoir can be switched off 
during the excitation process. In a real experiment this 
can be done, e. g., v ia optical or magnetic Feshbach res- 
onances [10, l4lTl42| , eventual residual finite interactions 
(e.g., due to magnetic field fluctuations for magnetic Fes- 
hbach resonances) can lead to a change of the excitation 
profile and thus represent a possible source of imperfec- 
tion. From numerical simulations, however, we find that 
for r <C J , the final temperatures (see Fig. [TUbD and 
the fraction of atoms in the dark state (see Fig. [TUb)) 
do not change significantly. From eq. (|33[) we find that 
T ~ 10~ 3 J 1 can be achieved for p^ = 5 x 10 14 cm~ 3 and 
a s = 20ah for 87 Rb in the lattice and 23 Na in the reser- 
voir. 



MANY PARTICLE COOLING 



In Fig.JSJj) and c) we plot the temperature of the lattice 
atom, as defined in eq. (|3"T|) and the height of the peak at 
q = against time in units of J° , where circles denote nu- 
merical results from the quantum trajectory simulations, 
and solid lines arc the analytical results from Levy statis- 
tics. In both cases we find excellent agreement of the nu- 
merical results with the analytical predictions. Typical 



In this section we will adapt the cooling scheme to 
many quantum degenerate bosons or fermions. We will 
assume that no interactions between the lattice atoms are 
present during the cooling process. In the case of bosons 
this means that the interaction Hamiltonian Hi = (see 
eq. p6p ). which can be achieved experimentally, e.g ., via 
an optical or magnetic Feshbach resonance ( [40l. |4U l42j| ) . 
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We consider a single species of Fcrmions, and thus s-wave 
interactions are forbidden due to Pauli blocking. 

The cooling process again follows the protocol as in the 
case of a single lattice atom: (i) In the excitation step 
Raman laser pulses are designed to excite atoms outside 
the dark state region to the first excited Bloch band, the 
coupling of the lattice atoms to the reservoir atoms are 
switched off during this step, (ii) In the decay step the 
dissipative coupling to the BEC reservoir randomizes the 
quasi-momcntum of lattice atoms decaying to the lowest 
band and consequently to a finite probability of falling 
into the dark state region. 

For non-interacting lattice atoms, the dynamics of the 
cooling process is again described by the master equation 
(|30[) . However, exact numerical simulations of the mas- 
ter equation are impractical, as the discretisation of the 
momentum space grid must be very fine to make pos- 
sible calculation of the low final temperatures. There- 
fore, we perform the analysis of the cooling process by 
numerical simulations of a quantum Boltzman master 
equation (QBME) 0, Ej|. The QBME is one of the 
simplest versions of the more general quantum kinetic 
master equation (QKME) [441 ] . which represents a fully 
quantum mechanical kinetic theory for the time evolu- 
tion of the system density operator, and was originally 
developed to analyse formation of Bosc-Einstcin conden- 
sates in atomic gases. The QBME is an equation for the 
diagonal elements w m = (m | p \ m) of the reduced sys- 
tem density operator, which describes the time evolution 
of the system in terms of classical configurations w m of 
atoms occupying momentum states m = [{m°} ? , 
in the two Bloch bands. Here denotes the occupa- 
tion of momentum state q in Bloch band a. In addition 
to the Born approximation and the Markov approxima- 
tion made in the derivation of the master equation, the 
QBME neglects the off-diagonal coherences between dif- 
ferent classical configurations contained in the QKME. 

In our numerical simulations we require a full QBME 
only for the decay step (ii), as for non-interacting atoms 
the excitation step (i) can be exactly calculated from the 
Heisenbcrg equations (|34p . i.e., from the excitation prob- 
ability Pj(q) as 

W tn\ = P j{l - fy/Mr*^ . ■ ( 38 ) 

The QBME for the second step (ii), the decay of excited 
lattice atoms back to the lowest Bloch band can be writ- 
ten as (see Appendix [C| 

w m =^T fc [m°_ fc (l±mj)«w -mj(l±m°_ fe )w m ] , 

k.q 

(39) 

where m' = m e-q-k,q and 

e q -k, q = [0, ...0, q i\ 0, ...0, -1, 0...0], (40) 
the upper (lower) signs are for bosons (fcrmions). 



Finally, we want to remark that the coherences which 
are neglected in the description of the decay step in terms 
of a QBME could, in principle, even be destroyed in a 
real experiment. For example, this could be done, by 
modulating the lattice depth and thereby randomizing 
the off diagonal elements, similar to the twirl in state 
purification protocols [46| . 

In the following we will present the results obtained 
from our numerical simulations first for the case of non- 
interacting bosons, then we will describe how finite in- 
teractions affect the excitation steps and finally we will 
present the results for spin polarized fermions. 

A. Results 

1. Non-interacting bosons 

For N non-interacting bosons, the T = ground state 
is the fully occupied q = momentum state, as in the 
case of a single atom. As a consequence we can use the 
same excitation pulse sequences as for a single atom (c.f. 
Fig. [H]). In Fig. \TTh ) we plot temperature against time 
which we obtain from numerical simulations [44j of the 
QBME cq. (f3"9"| . Temperature here is calculated by a 
Gaussian fit to the momentum distribution, excluding 
the central peak. Cooling to similar temperatures as 
in the single particle case can be obtained on shorter 
timescales for many atoms. This is due to the bosonic en- 
hancement factor, which appears as the factor (l + m°_ fe ) 
in the QBME (f3"9")) . In Fig. [TTb) we show the increase of 
the central peak against time. 




io 1 sxioVt 

FIG. 11: Numerical simulation of the QBME for N = 51 
bosons in M — 101 lattice sites, a) Temperature in units of 
4 J° against time in units of 1/ J° . b) Increase of the fraction 
of atoms in the dark state region qd £ [—0.06, 0.06] against 
time in units of 1/ J° . Parameters used as in Fig. [8] 



2. Excitation profile for bosons with interactions 

Until now we we have assumed that the interactions 
between the Bosons a is negligible. In this section we 
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investigate how small finite interactions will alter the ex- 
citation profile, and give approximate values for inter- 
action strengths that can be safely neglected. Here, we 
must compute the time evolution of the many-body sys- 
tem during the excitation step, which is described by the 
full two-band Bose Hubbard Hamiltonian for interacting 
atoms given in section UlI Al as Hq + Hls + Hi- 

This is possible in ID at temperature T = using time- 
dependent DMRG methods [43, H El- These numerical 
methods allow us to compute the time evolution of a ID 
many-body system where the Hilbert space can be ex- 
pressed as the product of local Hilbert spaces of dimen- 
sion d forming a ID chain. In these methods, the state 
of the system is written as a truncated Matrix Product 
state representation [47|, in which x states are retained 
in each Schmidt decomposition of the system. The dy- 
namical evolution for ID systems with sizes similar those 
seen in experiments can then be computed, starting both 
from the ground state and from weakly excited states. 

We simulate the time evolution in this way for the du- 
ration of one coherent Raman pulse. We consider the sit- 
uation with all the interactions equal, i.e. U Q0 = U 01 = 
U . To minimize the influence of box boundary condi- 
tions we used M = 41 sites and an initial Fock state of 
the form | ■ • • 0111110 • ■ - 0) with N = 5 atoms located 
at the centre of the system in the lower band. Since Fock 
states have a flat equally occupied momentum distribu- 
tion they allow the excitation profile to be determined by 
examining the final momentum distribution. A x = 50 
and an occupancy cut-off of n° lax < 4 and rc^ax — 2 
atoms per site for the lower and upper band respectively 
(equivalent to d = 15) was found to be sufficient. 

In Fig. [T2k ) the momentum distribution for the lowest 
band n°(q) after the pulse has been applied is displayed 
for a sequence of interaction strengths U aa /J°. For the 
parameters chosen when U aa = the initially flat mo- 
mentum distribution has a hole carved out at q = 4Aq, 
where Aq = lit/Aid. This demonstrates the selective ex- 
citation of this momentum and crucially the dark state 
at q = is left unchanged. As interactions are switched 
on the pulse begins to distort the final momentum distri- 
bution from this characteristic shape with a peak in the 
population of momenta either side of q = AAq as well as 
an increase in the population remaining at q = 4Aq it- 
self. This latter quantity provides a useful measure of the 
shape of the excitation profile and is plotted in Fig. WZb) 
and displays a linear increase only after U aa / J° > 0.2. 
The most relevant quantity for the cooling scheme is 
n (0), the population at q — 0, which is seen to increase 
linearly for small interaction strengths in Fig. 112b ). In 
total this reveals that for U aa <C 1/r <C | J x | the excita- 
tion profile remains close to that assumed for U aa = 0. 
If we include the additional constraint that U aa <C | J°| , 
which ensures that interactions do not substantially alter 
the ground state, the conclusions for the excitation pulse 
used earlier should not change. 
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FIG. 12: a) The momentum distribution in the lowest Bloch 
band obtained after a single excitation pulse applied is to an 
initial Fock state of N — 5 atoms in M — 41 lattice sites 
for a variety of interaction strengths U aa / J° € [0, 1.3]. The 
parameters used where Sqd = 1.23, w = 1.05J , (S — uS) = 
27.9J , V = lOWij, ujr = 2tt x 3.8kHz. b) The occupancy at 
momentum q — 4Aq and c) at q = against the interaction 
strength as a fraction of the flat occupation ng at = N/M. 



3. Fermions 

In contrast to bosons, in the case of N non-interacting 
(spin-polarized) fermions the T = ground state is char- 
acterized by the T = Fermi distribution, i.e., a step 
function n (q) = Q(q — <jf) for the momentum distri- 
bution in the lowest band, where the Fermi momentum 
qp = n(N — 1) /Aid. The excitation pulse sequence thus 
has to be changed in order to create a dark state region 
for atoms with momenta \q\ < qF- As a consequence 
the use of time square pulses is no longer advantageous, 
due to the large sidelobes they create in the excitation 
probability. We thus use a sequence of Blackmail pulses 
[3l| where these sidelobes are suppressed, as shown in 
Fig. 1131 where we compare typical excitation probabili- 
ties for a Blackman pulse and a time square pulse. The 
sequence of excitation pulses is now designed as in the 
previous cases, however the excitation probability can 
only be calculated numerically and the 7r-pulsc condition 
changes to r = 7r/0.42f2 for Blackman pulses [37j . 



In Fig. 114b ,) we plot the momentum distribution (again 
obtained from monte carlo simulations of the QBME) af- 
ter j = 0, 1, 2, 20 cooling cycles and find that the expected 
shape close to the expected T = Fermi distribution ap- 
pears after a few cooling cycles. Temperature is now 
obtained by fitting a Fermi distribution to occupation 
of the momentum states in the lowest Bloch band. In 
Fig. 114b ) we plot temperature against time in units of 
J° and find that temperatures UbT/AJ ~ 10~ 2 can be 
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FIG. 13: Comparison of the excitation probability for a 
time square (solid line) and a Blackman pulse (dashed line), 
shaded area depicts the filled Fermi sea. Parameters used: 
fi = 4.21J , 6qd = 0.84, (5 - u) = -27.6, M = 101, N = 51. 



obtained in tJ° ~ 500 for the parameters used (see cap- 
tion of Fig. [T4]) . The timescales for many fermions are 
slower than those of bosons due to Pauli blocking, which 
increasingly slows down the decay into an increasingly 
filled Fermi sea. 



n 



0.5 • 



a) 














r-r— 




n qd 



FIG. 14: Numerical simulation of the QBME for N = 
51 fermions in M — 101 lattice sites. a) Devel- 
opment of a sharp Fermi distribution after (dashed 
line), 1 (dash-dotted line), 5 (dotted line) and 20 (solid 
line) cooling cycles. b) Temperature in units of 4J° 
against time in units of 1/J°. Parameters used: fi = 
(11.58, 11.58, 1.1, 1.1) J°, {Sqd = 0.31,-0.31,1.2,-1.2), (5 - 
to) = -(28.19, 28.19, 26.89, 26.89) J° . 



VI. REALIZING COLD 
STRONGLY-INTERACTING GASES 

In this section we investigate how strongly correlated 
regimes can be realized where the cooling scheme can- 
not be applied directly This can be done by decoupling 
the optical lattice from the reservoir b and adiabatically 
ramping up the interaction strength. In this case our at- 
tention is restricted to the evolution governing by Hq+Hi 
with no Raman coupling and a = 0, and so describes the 
Bose-Hubbard model of the lowest band only. We anal- 
yse this by again making use of time-dependent DMRG 
methods |48|,|4j|. 



A. Ramping the ground state 

In principle adiabaticity requires an infinitely slow 
ramping. Here we determine a finite timescale in which 
near-adiabatic ramping can be achieved. We consider a 



ID lattice of M = 21 sites containing N = 10 atoms ini- 
tially prepared in the ground state with U 00 = and then 
raise the interaction strength over a time r r according to 



1 + exp 



t - r r /2 

T r /w 



(41) 

The constant C is fixed so as to make U 00 (t = 0) = and 
we choose the parameter w — 18 r r so that U 00 (t = r r ) ~ 
t^max- As this simulation begins from a non- interacting 
limit we took the occupancy cut-off to be n^ ax < 6 atoms 
per site (equivalent to d = 7). We found that for the 
computation of groundstates retaining x = 50 states in 
the matrix product decomposition was sufficient for this 
system, however, to accurately describe the dynamical 
ramping a \ = 100 was required. For the ramping itself 
we took J7 max = 20 J° and so approached the hard-core 
Bose lattice gas in ID (Tonks gas). 

To quantify the achievement of near-adiabatic ramp- 
ing we computed the energy deposited into the system 
£{T r ) — [E{t t ) — E ]/NJ° and the many-body overlap 
J-(r r ) = | ( ip(r r ) | -00 ) I as a function of the ramping time 
t t of the final ramped state | ip(r r )} with energy E(r r ) 
and the groundstate | tpo) with energy Eq for U 00 (r r ). 
Note that we express all energy differences £ as a frac- 
tion of NJ°. This is a useful energy scale since it is (up 
to a constant prefactor of order 1) the maximum energy 
difference which N non-interacting atoms inside the low- 
est band could have. Both £(r r ) and F(T r ) are plotted 
together in Fig. [T5k). We observe that for r r > 10/J° 
there is an exponential decrease in £(r r ) with r r , and 
thus for r r of order U 00 (T r )/J° there is negligible heat- 
ing within the Bloch band. Additionally F(r r ) can be 
seen to rapidly approach unity on the same timescale 
rigourously verifying that this final state is converging to 
the groundstate in the strongly correlated regime. We 
note, however, that since these calculations are based on 
a single system size they do not address the issue of scal- 
ing with M. Despite this we expect the results presented 
to be applicable to systems which are currently realized 
in experiments since their size is typically of the same 
order as considered here. 



B. Ramping weakly excited states 

In the previous section the degree of excitation induced 
by the ramping process on the ground state was quanti- 
fied. Here wc confirm that the near-adiabatic timescale 
determined for the ground state also applies to good ap- 
proximation to low- lying excited states. This is done 
by performing the ramping with a weakly excited initial 
state and demonstrating that the resulting final state re- 
mains weakly excited in the strongly interacting regime. 
Specifically, we generate an excited state by evolving the 
ryoo _ q g rounc j state for a time i ex in the presence of a 
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FIG. 15: a) The normalized energy difference £(r r ) (left-axis 
and dashed line) and the many-body overlap T(r r ) (right- 
axis and solid line) as a function of T r for a ramping of the 
f7 00 (t = 0) = groundstate with N = 10, M = 21 to 
U (t t ) = 20 J° according to the profile given in Eq. 04] 
b) The normalized energy difference £ (t) with the instanta- 
neous ground state at time t of the ramping with r r = 100/ J° 
for initial states generated with J° i ex equal to (i) 1, (ii) 5 and 
(iii) 9 respectively. These weakly excited SF states are char- 
acterized by £ = {0.01,0.02,0.03}, V = {1.2%, 5.7%, 7.6%} 
and T = {0.95,0.69,0.54} respectively for (i), (ii) and (iii). 
All three initial states had A£ ~ 0.06. 



spatially homogeneous on-site interaction which is vary- 
ing randomly in time in the range U 00 (t) £ [0, \ J }. With 
this method we constructed three excited initial states 
using </°t GX equal to (i) 1, (ii) 5 and (iii) 9. We charac- 
terize these states by their normalized energy difference 
£ , many-body overlap T with the ground state, quan- 
tum depletion T> (equal to N minus the largest eigen- 
value of the single-particle density matrix), and energy 
spread defined by (A£/NJ ) 2 = (H 2 ) - (H) 2 (see cap- 
tion of Fig. [T5] for specific numbers). Despite this being 
classical noise which coherently excites the ground state 
it does produce excited states with features similar to 
those found in experiments. In particular the charac- 
teristics for state (i) may be representative of a weakly 
excited SF state resulting from the cooling scheme out- 
lined. 

The three initial states were then ramped in an identi- 
cal way to section FVl Al using i> = 100/ J°. In Fig. 115b) 
the evolution of energy difference £ (t) as a function of 
the time t during a ramp is shown for each of them. 
From the two most excited of these states, (ii) and (iii), 
£ (t = t, ) is seen to level-off at around 4 times their 
initial values giving on the order of 10% heating within 
the Bloch band, and have A£(t = r r ) ~ 0.17, while the 
many-body overlaps of their final ramped states with the 
strongly-correlated ground state reduce to T = 0.55 and 
T = 0.41 respectively. This indicates that the ramping 
is not entirely adiabatic for these states. For the least 
excited initial state (i) we find that both £ (t = r r ) and 
A£{t = r r ) have approximately doubled, but crucially 
the heating £ (t = r r ) is still less than the 4%. Finally, 
for (i) the many-body overlap T = 0.89 suffers a smaller 
reduction and importantly remains sizable sizable show- 
ing that the final ramped state is still weakly excited in 
the strongly-correlated regime. 



VII. SPATIAL SIDEBAND COOLING IN A 
HARMONIC TRAP 

Up until now we have discussed the implementation of 
a Raman cooling scheme in an optical lattice with a flat 
external potential, based on the coupling of atoms in a 
lattice to a reservoir gas. As an additional remark, we 
describe in this section how these ideas can be applied 
in a different way to cool atoms in a harmonic trap to 
sites of lower energy, in analogy with sideband cooling to 
lower motional states in an ion trap. [351 ] . This example 
serves to strengthen the formal analogies between open 
quantum systems encountered in quantum optics, and 
systems of atoms in optical lattices coupled to a phonon 
bath. It can also be applied in the context of preparation 
of a quantum register with Fermions [39f ? ] , in this case 
compacting the cloud of atoms as they collect in the cen- 
tre of the trap, at the sites of lower energy (see Fig.lTfjk). 
Here we outline how single atoms can be cooled to sites 
of lower energy, and then briefly discuss the application 
to many fermions. 

We consider the same setup discussed in section|TJ with 
an identical Hamiltonian, Eq. (O, except that here we 
impose an additional external trapping potential (nor- 
mally a harmonic trap), He = Ej^Wj where hi is the 
number operator on site i. In this scheme, Raman cou- 
pling of atoms from the ground motional state in each 
lattice site to the excited motional state is not pulsed, 
but continuously switched on. The interaction with the 
external reservoir gas b, giving rise to the master equation 
([50)1 is also identical, and is also continuously switched 
on, giving rise to onsite decay. We assume here that the 
lattice depth is large, so that we are in a limit where col- 
lective processes (analogous to superradiance and subra- 
diancc in phonon emission) can be neglected (see section 

unci) . 



The basic concept of this process is shown in Fig. ITfjb . 
We choose the detuning of the Raman excitation so that 
it is below resonance with the excited motional level. If 
the atom subsequently tunnels (with amplitude J 1 ) to a 
neighbouring lattice site, and then undergoes a decay, it 
can be transferred to neighbouring lattice sites. (Note 
that tunnelling in the lowest band J° <C J 1 , and is fur- 
ther suppressed by the potential offset between the lat- 
tice sites, due to which tunnelling in the lowest band is 
no longer a resonant process.) Transfer to sites of lower 
energy is favoured, however, because of our choice of de- 
tuning for the Raman excitation. In Fig. 116b . dashed ar- 
rows indicate an effective two-step process including the 
Raman excitation and tunnelling to neighbouring sites. 
We note that this process can be made resonant for a site 
of lower energy, whilst far detuned from a site of higher 
energy. 

This is similar to the concept of sideband cooling in 
ion traps. There, the goal is to cool the motional state 
of an atom in a single trap, by coupling to states in 
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FIG. 16: (a) Techniques involving spontentous creation of ex- 
citions in the reservoir can also be used to cool Fermions in 
an external trapping potential to sites of lower energy, com- 
pacting them so that in the centre of the trap we observe unit 
occupancy, (b) Atoms are excited to the first excited motional 
level in each well, by a Raman transition (with two-photon 
Rabi frequency Q) detuned to energies below resonance with 
the excited state. Atoms can then tunnel to neighbouring sites 
with tunnelling rate J , and can decay back to the ground 
motional state by interaction with the external reservoir gas. 
The dashed arrows show the effective result of the two-step 
excitation and tunnelling process, and we see that transitions 
to sites of lower energy are favoured because these are closer 
to resonance, (c) Results for a single particle being cooled in 
this manner: Plot of the mean position of the atom in a sin- 
gle quantum trajectory as a function of time, beginning with 
an atom displaced 12 sites from the centre of a harmonic po- 
tential, ej = 0.25j 2 J 1 (solid line), and root-mean-square dis- 
placement of the atom from the centre of the trap, averaged 
over an initial uniform distribution between sites j = — 15 and 
j = 15 (dashed line). Parameters used: £1 = J , T = 2J 1 , 
S = -4J 1 



from monte-carlo wavefunction simulations. The solid 
line in the figure shows the mean position of the atom as 
a function of time for a single quantum trajectory, which 
illustrates the cooling and heating transitions as the atom 
moves closer to and further from the centre of the trap re- 
spectively. We also show the root-mean-square displace- 
ment from the centre of the trap [i.e., \J jhj)/M)] 

when we begin with a uniform distribution for the ini- 
tial position. The finite final temperature is determined 
by competition between heating and cooling steps, as is 
clearly illustrated by the transitions in the example tra- 
jectory. 

This method can be simply extended to fermions, 
where double occupation of the ground motional state in 
any lattice site is prevented by Pauli blocking. For this 
case it is possible to derive a QBME in analogy to that 
in sec. [Vj and to simulate these equations. Compaction 
of the fermions into the centre of the trap is observed, 
producing a quantum register with one fermion on each 
lattice site [51] . An additional step will, in general, be 
required to remove any remaining atoms from the upper 
motional level at the end of the cooling process. In the 
present form, this scheme is not well suited for bosons, as 
it would be difficult to prevent many atoms from collect- 
ing on a single site, even in the presence of interactions. 
However, the existence of this analogy is again a strong 
demonstration that ideas from quantum optics can be 
used in the context of this new type of open quantum 
system, an idea that has many possibilities for future ex- 
ploration. 



which the electronic state is excited, but the motional 
state is reduced by one quantum (this is called the red 
sideband). We can draw a figure analogous to Fig. [16b . 
so that in comparison with sideband cooling, we have re- 
placed excited and ground electronic states with excited 
and ground motional states in each well, and the coupling 
to lower motional states in sideband cooling is replaced 
by coupling to neighbouring sites of lower energy. In ion 
traps the goal is to cool an ion to the lowest motional 
state, whereas in this proposal, sideband cooling leads to 
a spatial redistribution of atoms towards the centre of 
the harmonic trap. 

In analogy with sideband cooling, the cooling steps 
in this scheme (coupling to sites of lower energy) will 
more strongly dominate heating steps (coupling to sites 
of higher energy) when the energy offset between the 
sites is made larger, and when the detuning S ~ — e, 
where e is the energy offset between sites. Note that in a 
harmonic trapping potential, where the potential differ- 
ence between neighbouring sites varies, the implementa- 
tion should involve a sweep in the detuning, making the 
cooling more efficient in different parts of the trap as a 
function of time. 

In Fig. 116b we plot example results for a single atom 
being cooled in a harmonic trapping potential, obtained 



VIII. CONCLUSIONS 

We have analysed a new cooling scheme based on ideas 
borrowed from sub-recoil Raman cooling schemes, but 
where these ideas are now placed in the context of a new 
form of open quantum system, where atoms in an opti- 
cal lattice are coupled to a BEC reservoir. In our case 
this setup provides laser assisted sympathetic cooling, in 
which the final temperature of atoms in the lattice is 
not limited by the temperature of the reservoir. This is 
motivated by the practical requirement of achieving low 
temperatures in an optical lattice, which is important for 
simulation of many strongly correlated quantum systems. 
Our scheme indeed provides a possibility to achieve low 
temperatures, and from our model we predict that tem- 
peratures a small fraction of the Bloch band width can 
be achieved in this way. 

Equally importantly, however, this work opens ques- 
tions and possibilities on the implementation of open 
quantum systems in this new context. Here we have 
demonstrated how ideas from quantum optics can be 
applied in a many-body system where we have strong 
control over many parameters, especially interactions be- 
tween atoms, and thus our system-reservoir coupling. 
The possibility exists to extend these ideas to different 
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types of dissipativc Hubbard models, and quantum reser- 
voir engineering more generally. These ideas could be 
used both in the context of state preparation and as- 
sist the engineering of new models with systems of many 
atoms in optical lattices. 

In the opposite direction, such open quantum systems 
also offer possibilities for investigating effects that are 
of fundamental interest in quantum optics, here in pa- 
rameter regimes that can be very different from regular 
quantum optics experiments. One example that is clear 
in the current context would be possibility to study the 
analogue of supcrradiant or subradiant decay, using the 
lattice and BEC parameters to engineer the wavelength 
of the excitations created in the BEC, and therefore the 
strength of these effects. 
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APPENDIX A: ADIABATIC ELIMINATION OF 
THE EXCITED STATE 

In this appendix we will show how the Raman coupling 
as given in eq. (|12p can be derived from eq. p0[) for large 
detuning A 3> u, , £1 of the two lasers from the excited 
state and for a one dimensional optical lattice. Starting 
point arc the Hcisenbcrg equations of motion for the field 
operators for the lattice atoms in ground and excited 
state, 

Mx,t) = -&&$ a (x,t) + iA$ B (x,t), (Al) 



where we have introduced 

n t (x) = n 1 { x ,t) + e- iSt n 2 (x,t). 

Formal integration of cqs. (|A1|) yields 



(A2) 



Mx,t) = -U dsn* s (x)e iA ^- s H a (x,s), (A3) 



where we have assumed tp e (x,0) = and f2t<o = 0. Ex- 
panding the field operators in terms of Bloch wave func- 
tions gives 

# a (x, s)«^ ^(xy^-^'-^t), (A4) 



and inserting this into eq. (|A3|) we find 

i nt(x,t) 



q . a 



id I 



■ 2e« - Lu a 



n* 2 (x,t) 



A + 2e« 



0J o 



A«(t), (A5) 



after integrating the resulting equation, and we have as- 
sumed Qi(t), £^2 (i) an d Aq(t) to be slowly varying on 
the timescale 1/A. 

As the detuning A of the two Raman lasers from the 
excited state is the largest frequency, i.e. A 3> w a , S, e", 
we can neglect all terms rotating with A, and thus 



i> e (x,t) 



€(x) 



fit 



2A 



2A 



A a At). 



(A6) 



Inserting into the Hcisenbcrg equations (|A1[) we can then 
read off the effective lattice Hamiltonian in an interaction 
picture with respect to Hq as 



^il 2 



4A 



E 

q,Q 

| £ i? Q , Q '(*z)(e i ^ + 

q,a,a.' 



,+S)t 



3 i( 2E ?+s,- 2 <+ w <»- w <.'- 5 )t^(^4 a! )t J 4 



q+k : a' • 



(A7) 



Here, we have assumed two running wave Raman lasers 
with relative momentum 5q and neglected the overlap of 
Wannier functions in different lattice sites. Transforming 
back to the Schrodinger picture with respect to the lattice 
and only taking into account the lowest two Bloch bands 
we obtain the Hamiltonian as given in eq. (|12[) . 



APPENDIX B: DERIVATION OF THE MASTER 
EQUATION 

In this section we will derive the master equation for 
the reduced system density operator p for the decay of 
lattice atoms from the first excited Bloch band back to 
the lowest band as given in eq. ([3H)l . As described above, 
the BEC can be treated as a (three dimensional) T — 
reservoir, and the Born-Markov master equation in the 
interaction picture is given by [43[ 



Kt) 



-Trr 



dt' 



H int {t),[H iat {t'),p(t)®p R ] 



(Bl) 



where pn is the density operator for the BEC reser- 
voir, and TrR {} denotes the trace over the reservoir 
states. In a standard way we change to the variable 
t = t — t' and extend the integration to r g [0, oo), 
assuming the correlation time in the BEC to be much 
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shorter than typical system time scales. Furthermore, 
we use J °° dte^ £_eo ^ = Tr5(e — e ) and 



APPENDIX C: DERIVATION OF THE 
QUANTUM BOLTZMANN MASTER EQUATION 



5k.k' 



0, 



(B2) 



for an effective T = reservoir. In the rotating wave 
approximation (i.e. neglecting terms rotating with fre- 
quencies u>) we then find 



- ~ 2 ft 



d 3 xu>i(x)u>o(x)e 



ikx 



\2ckpcl - 4 



C k C k/° - P c k c k 



(B3) 



where we have furthermore utilized J ,!,/ 1 ! <C u> and 
neglected the terms involving the band structure in the 
5 function. Note that we have also neglected the energy 
shift terms resulting from principal value integrals in this 
derivation. The energy shifts resulting from two atoms 
immersed in a 3D reservoir have been computed to be 
very small [H^ . 

In order to specialize to a one dimensional optical lat- 
tice in the three dimensional reservoir we convert the sum 
in eq. (|B3[) into an integral in spherical coordinates, and 
integrate over |k| (which is fixed by energy conservation) 
and the azimuthal angle. For a one dimensional lattice 
we then obtain 



J2 -f (24pcj; 



c k c k p - pc k c k 



(B4) 



after converting the integral over the polar angle back to 
a sum over the momentum in the x-dircction, where Tk 
is given by eq. (|32[) in harmonic oscillator approximation. 

We can also rewrite the master equation in terms of 
position space operators, af = (l/VTI) J2 q exp(iqxi)A q . 
Defining ci = (aj-^aj We obtain 



r_ fc 

2 

fe 



+ E ~r E ( 25 ^ 5 I - 4 & ip - p4c* 



(B5) 



When kd ^> 27r, the terms in the sum over i ^ j decay 
very rapidly with \i — j\, as is shown in [39[ for the same 
master equation with a fcrmionic reservoir. In this limit, 
collective effects of decay involving atoms on different lat- 
tice sites can be neglected (i.e., there are no superradiant 
or subradiant effects). See the discussion in section UlI CI 
for more details. 



In this section we will derive the quantum Boltzmann 
master equation (QBME), eq. ([55]) which we use to ana- 
lyze the dissipative coupling of lattice atoms to the BEC 
reservoir in the case of many non-interacting bosons and 
fermions, as presented in section [Vj The QBME is an 
equation for the diagonal elements w m — (m | p | m) of 
the system density operator p and can be derived from 
the master equation ([50)1 . We project on the diagonal 
elements of the density operator and find 



k n, q,q' 



2<m| (A°_ fe )Ul| n) ( n \(A^A° q ,_ k \m) 
-<m|K)U° (4> fc M,|n) (n|m)« 



(CI) 
(C2) 
(C3) 



The expectation values in the first part (|C1[) of this equa- 
tion can be calculated as 



(m|n) (n|(Ai)U° fe (A° )UJ, |m) 



(n|(^)t^,_ fc |m 



m,n+e 



q'-k.q' •> 



(m| {A q _^A\ |n) = Jm q _ k Jl ± mj<5 m , n+ 



(C4) 



where the upper (lower) signs are for bosons (fermions), 
the product of the two S functions appearing in these ex- 
pressions can be evaluated as 5 m ,n+e q , _ k (J /<5m,n+e,_ fc , = 
^m.n+e,_ fc &q,q< ■ For the second and third line (|C2[) and 



(|C3[) we find the identical expressions 

m\(l ± m q _ k )8 q . q >5 m ^. (C5) 

Thus, we can perform the sums over q' and n and end 
up with the QBME 



E 

k . q 



r fc [ m °^fe( 1 ± mj)w m ' - mj(l ± m q _ k )w m ] 



(C6) 



with m' = m e„_ 



q-k,q- 
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